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Among  anti-submarine  warfare  (ASW)  target  motion  analysis 
(TMA)  algorithms  for  mobile  platforms,  passive  bearings  only  track- 
ing is  the  most  basic  segment.  The  present  state-of-the-art  of 
tracking  by  this  method  is  not  satisfactory  [ 1 ] . None  of  the  existing 
extended  Kalman  filter  algorithms  yield  both  accurate  tracks  and  an 
accurate  solution  quality  indicator.  Generally,  these  algorithms 
perform  no  better  than  much  simpler  manually  controlled  algorithms. 
Several  approaches  to  the  use  of  the  recursive  filtering  theory  are 
discussed  in  references  [l,  2].  The  approaches  can  all  be  refer- 
enced to  a basic  text  book  Extended  Kalman  Filters.  Differences 
among  the  algorithms  correspond  to  choice  of  coordinate  systems, 
whether  one  aligns  the  linearization  point  of  the  nonlinear  measurement 
to  the  line  of  sight  or  not,  whether  or  not  one  limits  artificial  range 
collapse,  and  whether  or  not  one  uses  some  other  form  of  divergence 
control.  The  difficulties  are  quite  well  documented  in  reference  [l]. 

The  severity  of  the  undesired  behavior  differs  from  algorithm 
to  algorithm.  This  difference  depends  on  the  nature  of  the  algorithm 
and  the  accuracy  of  the  measurements,  the  rate  of  the  measurements, 
the  range  and  data  rate,  and  the  geometric  configuration  of  OWN  ship 
and  target  ship.  However,  the  fact  that  none  of  the  algorithms  generally 
give  satisfactory  performance  leaves  one  to  question  whether  or  not 
the  difficulties  lie  with  the  nature  of  the  extended  Kalman  filter 
which  is  being  applied  to  this  problem.  It  is  the  purpose  of  this  paper 

=*This  work  was  partially  supported  by  the  Air  Force  Office  of 
Scientific  Research  under  contract  AF-AF0SR-44620-75-C-0023 
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to  point  out  the  difference  in  performance  between  a simple  extended 
Kalman  filter  and  an  optimal  Bayesian  filter  in  pictorial  form.  And, 
to  show  reasons  why  extended  Kalman  filters  do  not  perforin  satis- 
factorily. 

The  actual  implementation  of  the  Bayesian  optimal  filter  is  a 
Gaussian  sum  approach  documented  in  references  [2,  4,  5].  How- 
ever, the  purpose  of  this  paper  is  not  to  demonstrate  the  utility  of 
the  Gaussian  sum  approach  in  nonlinear  filtering,  but  instead  to  show 
the  difference  between  an  optimal  Bayesian  filter  implemented  in  any 
manner  whatsoever  and  a linearized  or  expended  Kalman  filter.  It  is 
hoped  that  the  following  description  will  give  the  reader  insight  into 
the  nature  of  both  the  optimal  Bayesian  filter  and  the  approximation 
nature  inherent  in  any  extended  Kalman  filter. 


2. 


PASSIVE  BEARINGS  ONLY  TRACKING  EXAMPLE 


A simple  two-dimensional  version  of  the  bearings  only  track- 
ing problem  is  considered.  The  state  vector  is  chosen  to  be  two- 
dimensional  (latitude,  longitude)  for  simplicity.  The  extention  to 
the  case  of  four  state  variables,  including  bearings  and  speed  in- 
formation is  straight  forward.  Again,  for  simplicity,  OWN  ship 
is  assumed  to  be  traveling  at  a constant  speed  around  a circular 
path.  Any  trajectory  for  OWN  ships  motion  could  have  been  used. 

These  simplifications  have  been  introduced  because  the  purpose 
of  this  paper  is  to  illustrate  the  nature  of  the  optimal  Bayesian  approach 
to  this  problem  and  the  effects  of  various  linearization  approximations, 
not  to  solve  the  problem  in  the  most  expeditious  manner. 

Once  an  engineer  sees  the  nature  of  an  "optimal"  solution  to  a 
nonlinear  filtering  problem,  he  can  often  develop  a simpler,  but 
adequate  suboptimal  (often  a linearization  type)  filter. 

The  state  vector  propagates  according  to  the  linear  plant  model. 

±k+l  + 1 

where  the  two-dimensional  state  represents  latitude  and  longitude 
coordinates. 


and  the  state  is  observed  by  the  scalar  nonlinear  measurement 
function: 


Zk  = hk(-k)  + Vk 
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The  above  model  arises  in  connection  with  the  tracking  geometry  of 

T 

Figure  1,  where  the  target  T at  the  position  defined  as  = (x^,  y^) 
is  undergoing  a random  walk  in  the  two  dimensional  state  space.  The 
observer  S is  passively  measuring  the  line  of  sight  or  it  travels  in 


Figure  1.  Geometrical  Definition  of  Vector  Tracking  Example 


a deterministic  orbit  around  the  unit  circle.  It  is  important  to  note 
that  when  P is  equal  to  zero  the  target  is  unobservable  and  the 
filtering  problem  becomes  degenerate. 

The  problem  to  be  considered  here  is  how  to  obtain  the  best 
possible  estimate  of  the  state  x,  of  the  system  (1-2)  based  on  all 

k —k 

data  up  to  time  k,  z . * 


*The  collection  of  data  up  to  time  k jf  is  defined  as 

*k. 
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3. 


BAYESIAN  RECURSION  REIATIONS 


Since  the  state  vector  x,  is  a random  variable,  in  the  prob- 
abilistic  context  of  this  discussion,  the  a posteriori  density  p(x^|z  ) 
provides  the  most  complete  description  possible  of  x^.  This 
density  is  determined  recursively  from  the  Bayesian  recursion 
relations  given  below.  The  filtering  density  is  given  first  and  the 
prediction  density  is  next. 


p(*kUk) = ckp{xh^zk"1)p(zk^xk) ' 

where  the  normalizing  constant  c is  given  by 

K 

1 /ck  = p(zk,zk_1)  =XnP(X k,2kl)p(zk,xk)dxk’ 


and  the  initial  condition  for  (8)  is 


p(x  | z ) Ap(x  ) = N(y  -x'p') 
o = o o o o 


The  densities  p(z^|x^)  and  p(x^|x^  j)  are  determined  from  (1) 
and  (2),  and  the  a priori  densities  p(v^)  an(3  p(w 

p,*kK-l)  = N|xk'xk-rQk) 


A 1 | A t -1  * | 

:-x,A)A  : io"rt/-T  exp  < 5(x-x)  A (x-x)  > 

|a|  I/2u^/z  1 ’ 


4. 


THE  EXTENDED  KALMAN  FITTER 


The  nonlinear  filtering  problem  is  solved  when  the  density 
p'kv  can  be  obtained  for  all  k . However,  except  when  equation 
2 is  linear  and  the  a priori  distributions  are  Gaussian,  it  is 
generally  impossible  to  determine  p(x^|  z^)  in  a closed  form  using 
(8)  — (1 1).  In  the  linear  Gaussian  case,  the  relations  describing  the  condi- 
ditional  mean  and  covariance  of  p(x^|  z^)  are  the  well-known  Kalman 
filter  equations  (6).  The  difficulties  associated  with  the  explicit 
determination  of  the  a posteriori  density  have  led  to  the  development 
of  approximate  procedures  for  estimating  the  state  of  nonlinear 
stochastic  systems.  The  most  commonly  used  approximation 
involves  the  assumptions  that  the  a priori  distributions  are  Gaussian 
and  that  the  nonlinear  system  can  be  linearized  relative  to  latest 
estimate  yielding  the  "extended"  Kalman  filter  which  has  seen  wide- 
spread application  to  nonlinear  systems. 

The  extended  Kalman  filter  has  performed  well  in  many 
applications  but  there  are  numerous  examples  in  which  unsatis- 
factory results  have  been  obtained.  This  has  spurred  the  develop- 
ment of  other  procedures.  Most  of  these  either  implicitly  or 
explicitly,  retain  the  assumption  that  p(3c^ I z^)  is  Gaussian  and 
essentially  provide  a means  for  modifying  the  mean  and  covariance 
of  the  density.  However,  the  Gaussian  assumption  greatly  reduces 
the  generality  or  richness  of  possible  a posteriori  densities  by  throwing 
out  all  possibility  of  more  realistic  multimodal  densities. 

Starting  at  stage  0 and  linearizing  about  a point  in  state 

space  and  dropping  all  but  the  first  order  terms,  we  can  solve  the 
Bayesian  recursion  relations  to  find: 


°k  = H«ik»  k HT,i> + \ 
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r » 


i 


\ ■ pL  H’r,v a:' 

15 

\ ^ + Kk|zk  - h&» 

16 

pk  = - Kk  pL 

17 

pUi  = pk + °k 

18 

a / * 

19 

The  various  versions  of  the  extended  Kalman  iilter  are  simply  dif- 
ferent choices  of  x^.  The  most  common  text  book  choice  is 

Xk' 

H(xr)  = &h(x)/bx|^_  . 


Iterated  extended  Kalman  filters  require  processing  the  data 
through  equations  (14),  (15),  and  (16)  two  or  more  times  while  only  up- 
dating x^  > starting  with  x^  equal  to  x^  . If  the  filter  performs 
in  a reasonable  manner  this  iteration  will  move  x^  toward  a point 
x*  such  that 


Zk  = 


20 


which  is  "near"  the  original  point  x^  in  some  sense.  The  difficulty 
with  linearization  about  a point  that  does  not  satisfy  equation  (20)  can 
be  seen  from  the  following  simple  example. 

Consider  a case  where  the  predicted  state  vector  x^  is 


i-C) 


J 


8 


rind  the  measurement  equation  is  of  the  form: 


zk  = W + vk  = ' 1 *k  + |xk-yk'  + 
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I I 

' 


p(vk)  = N(vk>  . 1) 


For  this  measurement  function: 
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p,zk>V  yk)  = N(zk‘-lxk-(xk-yk»-  •*> 


23 


Using  the  linearization  procedure  of  the  extended  Kalman  filter  and 

linearizing  about  x^,  i.e.,  (0,  0)  we  have: 

3 h^  (0,  0)  2hk(0,  0) 


or 


\ < v V = \ (0’  °>  + — + 


\ <v  V = °- 1 yk 
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and  using  this  to  replace  h^  in  the  measurement  density  the  density 
used  by  the  Kalman  filter  to  represent  equation  (23)  is: 


PKAl.‘Wyk>*  Nl\--1Xk'  •» 
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The  true  function  and  this  usual  approximation  to  it  are  shown  in 
Figures  2a  and  2b,  respectively  for  the  case  when  7^  is  taken  equal 
to  one.  The  fact  that  the  nonlinearity  is  quite  large  with  respect  to 
the  linear 


2 Gaussian  Approximation 


term  is  born  out  by  the  fact  that  the  linearization  about  the  value 
of  states  predicted  before  the  measurement  leads  to  a very 

poor  measurement  density  function  approximation  and  thus  to  a 
very  poor  a posteriori  density  function  approximation. 

In  many  cases  it  is  preferable  if  the  function  h^  can 
be  linearized  about  a point  when: 


zk'W  yk)  = 0 26 

This  region  is  described  by  the  equation: 

\ - ■ lxk  ' ,xk  - yk>2  = 0 17 

with  z,  = 1 . 

k 

Which  is  essentially  a U shaped  region  extending  to  infinity  in  two 
directions.  This  region  is  shown  for  a limited  space  in  Figure  2a  as 
the  U shaped  ridge  on  the  true  measuiement  density. 

Consider  the  iteration  procedure  which  picks  out  values 
closest  to  the  previous  expected  value  or  mean  of  the  state  and 

satisfying  (27).  Note  that  all  values  of  and  y^_  such  that: 


W = yk  i 

- • lxk 

28 

W c yk  - 

. 05  + Jzk  + . 0025  - . 05  yk 

29 

satisfy  equation  (27).  Using  this  and  replacing  the  value  appearing  on 
the  right  hand  side  of  these  eqiiations  by  the  previous  mean  four 
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0 

Figure  3.  Points  Satisfying  Equation  (27) 

In  linearizing  the  nonlinearity  about  each  of  these  four  points  the 
four  resulting  functions  reduce  to  two  giving  the  two  approximation 
to  p(z^  !x^)  shown  in  Figure  2c.  It  is  seen  that  this  is  a considerably 
better  approximation  to  the  true  function  than  that  of  Figure  2b,  and 
that  it  might  be  improved  still  more  by  an  increase  of  the  predicted 
variance  of  each  of  the  Gaussians.  It  should  also  be  noted  that 
iterating  the  linearization  to  its  limit  point  will  lead  to  the  right 
or  left  hand  ridge  of  Figure  2c,  which  is  considerably  better  than 
the  non-iterated  case  but  which  ignores  the  second  ridge  completely. 

An  alternate  approach  would  be  to  require  that  the  linearization 
occur  at  the  point  nearest  xf  satisfying  equation  (26).  This  would 

K 

only  be  acceptable  if  the  measurement  data  was  perceived  to  be  more 
accurate  than  previous  state  information. 
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5. 


GAUSSIAN  SUM  APPROACH 


When  the  linearization  approximations  discussed  in  the  last 
section  are  not  satisfactory,  one  must  return  to  the  general  Bayesian 
recursion  relations  of  section  3.  It  is  generally  not  possible  to 
solve  these  equations  analytically  and  some  numerical  or  approxi- 
mation technique  must  be  used.  These  include  the  point  mass  method 
of  Bucy  and  Senne  [7],  spline  approximations  by  Jan  and  Figueirado 
[8]  and  the  Gaussian  sum  methods  by  Alspach  and  Sorenson  [2,  5].  A 
good  summary  of  these  techniques  is  given  in  reference  [9]. 

The  Gaussian  sum  implementation  is  used  here.  However,  the. 
example  is  considered  more  important  than  the  methods,  so  that 
only  a brief  outline  of  the  methods  is  given.  For  a more  complete 
description  of  the  approach  the  reader  is  referred  to  the  references. 

The  ba.sic  approach  for  this  problem  is  to  approximate  the 
measurement  density  p(z^|x^)  by  a weighted  sum  of  Gaussian  or  a 
Gaussian  sum. 


p(zk^Xk>  = N(fk  ‘ h(^)’ 

M 

PA<zkK»=S“i  N<\-HkiV  V 

i = l 
N 

f a.  = 1 


a.  £ 0 


30 


31 


32 


i=  1 


f i ;he  example  in  section  4,  the  two  ridges  in  Figure  2c  could  be  a 
tv  j term  Gaussian  sum  approximation  to  the  function  of  Figure  2a. 

A more  detailed  approximation  would  be  that  of  Figure  2d  where  a 
30  term  Gaussian  sum  approximation  to  the  function  of  Figure  2a  is 
used. 
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t.V;  a" 


A 


Given  a Gaussian  or  a Gaussian  sum  expression  for  p(x^  | z ) 
in  equation  (8)  and  the  Gaussian  sum  approximation  to  p (z^  !x^  ) from 
equation  (31),  one  obtains  an  analytical  Gaussian  sum  expression  for 

Pf^Jz  ) from  equation  (8).  Equations  (9)  and  (10)  can  also  be  solved 

i k 

analytically  given  a Gaussian  sum  expression  for  pfx^j  I z )•  This 
can  be  repeated  to  give  a complete  solution  for  the  a posteriori 
density  function  based  on  any  measurement  realization. 

I he  only  reason  this  approach  does  not  lead  to  a true  repro- 
ducing density  is  that  the  number  of  terms  in  the  Gaussian  sum  for 
pfx^lz  ) is  larger  than  the  number  in  p(x^  Jz  ” ) if  no  other 
approximations  are  introduced.  This  is  discussed  at  length  in 
the  reference  where  two  additional  approximations  are  introduced 
and  analyzed.  The  first  is  that  if  any  term  in  the  Gaussian  sum  has 
a weighting  o of  less  than  £ it  is  dropped.  The  second  approxi- 
mation is  that  if  two  terms  in  the  sum  become  similar  in  an  sense 
they  are  combined  and  their  weighting  coefficients  (a  , a ) are 

, . 1 j 

combined. 


i.  NUMERICAL  RESULTS 

The  geometry  of  the  situation  is  repeated  in  Figure  4a.  Then 
a Gaussian  sum  approximation  is  found  to  approximate  the  measurement 
ment  density  function  p ( J ) which  can  be  written  as: 

z z 

= const  exp  [-.  5(2^  - h^fx^))  /avJ 

The  true  measurement  function  is  shown  in  Figure  4b  where  it  has 
been  assumed  that  the  measurement  density  also  contains  the  infor- 
mation that  the  probability  is  zero  that  the  target  is  greater  than  6 
orbital  radii  away  from  the  observer.  This  accounts  for  the  sharp 
cutoff  seen  in  the  figure.  Eesidcs  this  cutoff  the  major  feature  of 
this  function  that  distinguishes  it  from  any  possible  one  Gaussian 
approximation  is  the  cone  shape.  This  spreads  out  away  from  the 
observer  showing  that  the  farther  the  target  is  from  the  observer 
the  larger  the  possible  absolute  error  in  the  targets  position,  h or 
this  example  the  additive  measurement  noise  has  a one  sigma  value 
of  . 1 radians  or  about  6 degrees. 

'When  the  basic  extended  Kalman  procedure  is  used  this 
function,  (h^),  is  linearized  about  the  previous  best  estimate  for 
state.  "When  this  linearization  is  performed  about  any  point  away 
from  the  line  of  sight  of  the  latest  measurement  (z^  - 1 in  this  case) 
very  bad  results  can  occur  as  showm  in  the  example  of  section  4. 

These  equations  were  implemented  by  Bucy  and  Senne  for  this  example 
in  reference  [l]  and  divergence  occurred.  In  Figure  4c  the  non- 
linearity has  been  linearized  about  a point  on  the  line  of  sight  of 
the  measurement  and  at  a distance  of  3 units  from  the  observer. 

This  gives  the  correct  value  for  the  variance  at  the  point  about  which 
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Targct/OWN  Ship  Tracking  Geometry 


the  linearization  is  made,  too  email  a value  for  pointe  farther  from 
the  observer,  and  too  large  a value  for  pointe  closer  to  the  observer. 
It  can  also  be  seen  from  Figure  4c  that  there  is  no  way  to  incor- 
porate the  cutoff  data  in  the  one  Gaussian  approximation  in  a simple 
fashion.  It  is  seen  that  this  density  function  runs  to  infinity  in  both 
directions  in  the  state  space  plane.  Figure  4d  shows  a 10  Gaussian 
sum  approximation  to  p(z^|x^)  of  Figure  4a  which  is  obtained  by  the 
methods  of  reference  [2,4]  and  in  which  no  search  was  required. 


This  technique  is  applied  to  a dynamic  example  in  Figure  5 
where  the  position  of  the  observer  is  shown  by  the  cross  on  the  unit 
orbit  and  the  cross  on  the  density  function  shows  the  true  position  of 
the  target.  The  a priori  estimate  for  the  initial  state  was  taken  to 
be; 


while  the  true  value  ofthe  initial  state  (and  all  subsequent  values 
since  there  is  no  plant  noise)  was  taken  to  be; 


The  value  of  p on  the  curve  is  the  value  of  p (k- 1 ) where  p is  10 
degrees  and  p^  is  -90  degrees.  The  measurement  noise  here  has  a 
one  sigma  value  of  . 01  radian  or  about  one-half  degree.  The  non- 
gaussian  a posteriori  filtering  density  function  is  seen  to  propagate 
from  stage  1 to  stage  9 in  this  figure  where  a measurement  is  taken 
every  10  degrees.  There  were  6 terms  in  each  Gaussian  sum  approxi 
rnation  to  p( | ) in  thiB  example  and  the  maximum  number  of 
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tirms  in  any  filierinp  density  was  117.  The  combination  and  droppinp 

-4 

criteria  were  chosen  to  be  f>  ^ = 6^  = 10 

Figure  6 shows  the  error  in  the  estimate  of  the  state  from 
the  Gaussian  s\im  filter  of  Figure  5 and  the  extended  Kalman  filter. 

In  addition,  a divergence  measuring  parameter  which  should  have  an 
average  value  of  2.  in  a Monte  Carlo  experiment  is  shown. 

In  Figure  7 the  measurement  one  sigma  value  Vt'as 
increased  to  . 02  radians  or  about  one  degree  and  a plant  noise  with 
one  sigma  value  of  . 1 was  added.  The  statistics  of  the  initial 
estimate  were  changed  to: 


Here  the  propagation  of  the  distinctly  non  Gaussian  filtering  density 
from  stages  1 to  6 is  shown. 

While  these  results  show  the  ability  of  the  Gaussian  sum 
approximation  to  calculate  non  Gaussian  a posteriori  densities  in 
a vector  case,  complete  verification  of  these  results  have  been 
made  by  comparison  with  the  considerably  more  expensive  results 
of  Bucy  [7]  and  by  a Monte  Carlo  simulation.  Plant  noise  was  added 
to  the  system  and  the  initial  conditions  and  a priori  statistics  were 
changed  to  be  consistent  with  the  Monte  Carlo  results  presented  in 
[7].  The  plant  and  measurement  noise  were  again  white  Gaussian 
sequences  with  O = 0.  1,  V#°-  1 rad/stage,  = I,  and 
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A Monte  Carlo  average  of  100  runs  was  performed  filtering  the 

sample  paths  for  the  Gaussian  sum  and  extended  Kalman  filters.  The 

results  are  presented  in  Table  1.  With  the  plant  noise  added,  the 

increased  measurement  noise,  and  the  higher  rate  of  rotation  of  the 

observer  in  its  orbit,  the  extended  Kalman  filter  performance  with 

respect  to  the  Gaussian  sum  filter  was  greatly  improved.  Note  that 
o 

this  was  not  true  of  the  more  basic  extended  Kalman  filler  imple- 
mented by  Bucy  and  Senne  [7],  which  was  linearized  about  the  latest 
estimate  for  state  and  not  about  the  latest  measurement.  That 
version  still  had  severe  divergence  characteristics.  Table  1 shows 
the  results  of  the  filtered  estimate  for  10  stages.  The  average  error, 
average  covariance,  and  average  divergence  parameter  are  presented. 
It  should  be  noted  from  the  average  covariance  that  the  optimal 
filter  continues  to  give  superior  performance  even  at  the  later 
stages.  the  divergence  parameter,  which  should  be  such  that 
A = E(A)  = 2 for  an  infinite  number  of  runs,  is  consistently  larger 
for  the  extended  Kalman  filter.  The  large  values  exhibited  by  this 
parameter  for  some  stages  result  from  sample  paths  for  which  the 
covariance  of  the  extended  Kalman  filter  becomes  nearly  singular.  It 
is  worth  nofmg  that  the  results  reported  here  correspond  closely 
with  results  obtained  by  Bucy  and  Senne  using  the  method  described 
in  [7J.  This  lias  been  established  in  private  discussions  in  which 
both  methods  were  applied  to  the  same  measurement  realizations  and 
with  all  other  assumptions  identical. 
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